import numpy as np
import struct
from matplotlib import pyplot as plt
import time

def sigmod(x):
    result = 1.0/(1.0 + np.exp(-x))
    return result


def softmax(x):
    x_exp = np.exp(x)
    # 如果是列向量，则axis=0
    x_sum = np.sum(x_exp, axis=1, keepdims=True)
    result = x_exp / x_sum
    return result

def DS(x):
    f = sigmod(x)
    result = f * (1.0 - f)
    return result

def DSoft(x):                               #这个对吗？
    #x.shape=(1,10)
    x_exp = np.exp(x)
    x_sum = np.sum(x_exp, axis=1, keepdims=True)
    #result = (x_sum - x)/x_sum**2
    Y = softmax(x)
    #Y.shape=(1,10)
    #print(Y)
    #Y是[[...]],需转换为[],故取Y[0]
    #Y' = np.diag(Y[0])
    #print(Y)#1，
    result = np.diag(Y[0]) - Y.T @ Y
    return result

def ReLU(x):
    if(x<0.0):
        reslut = 0.0
    else:
        result = x
    return result

def forward(x,w1,w2,b1,b2):             #正向传播

    a1 = x
    a2 = np.dot(a1, w1.T) + b1
    z2 = sigmod(a2)
    a3 = np.dot(z2, w2.T) + b2
    z3 = softmax(a3)

    return z2, z3, a2, a3
    
def loss(y_hat, y):
    #y.shape[0] 是神经元的个数
    result = np.array((y-y_hat)**2) / 2 #* y.shape[0]    为什么乘最后一层神经元的个数？
    result = np.sum(result, axis=1) / 10
    return result

def Dloss(y_hat, y):
    result = np.array(y-y_hat)             #只有y是变量，y_hat不是变量，对y求导，因此是这个
    #result = np.array(y_hat-y)            #为什么不是这个？
    #result = np.array(abs(y_hat-y))       #为什么不是这个？
    return result

def grad(out_hat, out, a2, a3, w1, w2, b1, b2, z2, x, xite):
    '''
    dw1                 mid*in
    dw2                out*mid

    Dloss(out_hat, out)   1*10
    DSoft(a3)            10*10
    w2                  10*mid
    DS(a2)               1*mid
    '''
    #实现参数更新
    #dw1 = (Dloss(out_hat, out) * DS(a3) @ w2 * DS(a2)).T @ x            #仔细算一下
    dw1 = (Dloss(out_hat, out) @ DSoft(a3) @ w2 * DS(a2)).T @ x
    #dw2 = (Dloss(out_hat, out) * DS(a3)).T @ z2
    dw2 = (Dloss(out_hat, out) @ DSoft(a3)).T @ z2

    #db1 = Dloss(out_hat, out) * DS(a3) @ w2 * DS(a2)
    db1 = Dloss(out_hat, out) @ DSoft(a3) @ w2 * DS(a2)
    #db2 = Dloss(out_hat, out) * DS(a3) 
    db2 = Dloss(out_hat, out) @ DSoft(a3) 

    w2 = w2 - xite * dw2
    w1 = w1 - xite * dw1
    b2 = b2 - xite * db2
    b1 = b1 - xite * db1

    return w1, w2, b1, b2


def decode_idx3_ubyte(idx3_ubyte_file):#此函数用来解析idx3文件，idx3_ubyte_filec指定图像文件路径
    #读取二进制数据
    bin_data=open(idx3_ubyte_file,'rb').read()
    #解析文件头信息，依次为魔数、图片数量、每张图片高、每张图片宽
    offest=0
    fmt_header='>iiii'    
    magic_number,num_images,num_rows,num_cols=struct.unpack_from(fmt_header,bin_data,offest)
    print('魔数：%d,图片数量：%d，图片大小：%d%d' % (magic_number,num_images,num_rows,num_cols))
    #解析数据集
    image_size=num_rows*num_cols
    offest += struct.calcsize(fmt_header)
    fmt_image='>'+str(image_size)+'B'
    images=np.empty((num_images,num_rows,num_cols))
    for i in range(num_images):   
        #if (i+1)%10000==0:
            #print('已解析%d'%(i+1)+'张')        
        images[i]=np.array(struct.unpack_from(fmt_image,bin_data,offest)).reshape((num_rows,num_cols))
        offest+=struct.calcsize(fmt_image)
    return images
'''images是一个三维数组,images[i][a][b]表示第i张图片的倒数第a行，b列的像素'''

def decode_idx1_ubyte(idx1_ubyte_file):#解析idx1文件函数，idx1_ubyte_file指定标签文件路径
    #读取二进制数据
    bin_data=open(idx1_ubyte_file,'rb').read()
    #解析文件头信息，依次为魔数和标签数
    offest=0
    fmt_header='>ii'
    magic_number,num_images=struct.unpack_from(fmt_header,bin_data,offest)
    print('魔数：%d，图片数量：%d张' % (magic_number,num_images))
    #解析数据集
    offest+=struct.calcsize(fmt_header)
    fmt_image='>B'
    labels=np.empty(num_images)
    for i in range(num_images):
        #if (i+1)%10000==0:
            #print('已解析：%d'%(i+1)+'张')
        labels[i]=struct.unpack_from(fmt_image,bin_data,offest)[0]
        offest+=struct.calcsize(fmt_image)
    #print(labels[0])
    return labels
'''labels是一个一维数组，每个元素都一一对应images[i]'''

def label_reshape(output):
    #output是原始的标签，内容为5，8等；out是更改为10个神经元的标签，内容为0,0,1,0,0,0,0,0,0,0
    output = output.reshape(-1,1)
    out = np.zeros((output.shape[0],10))
    for i in range(output.shape[0]):
        if output[i] == 0:
            out[i][0] = 1
        elif output[i] == 1:
            out[i][1] = 1
        elif output[i] == 2:
            out[i][2] = 1
        elif output[i] == 3:
            out[i][3] = 1
        elif output[i] == 4:
            out[i][4] = 1
        elif output[i] == 5:
            out[i][5] = 1
        elif output[i] == 6:
            out[i][6] = 1
        elif output[i] == 7:
            out[i][7] = 1
        elif output[i] == 8:
            out[i][8] = 1
        elif output[i] == 9:
            out[i][9] = 1
    return out, output


#5.45.2
#8.23
'''------------------------------------------主函数开始---------------------------------------------'''

#读取数据集
train_image = r"D:\VScode\project\python\神经网络实验\MNIST\train-images.idx3-ubyte"
train_label = r"D:\VScode\project\python\神经网络实验\MNIST\train-labels.idx1-ubyte"
test_image = r"D:\VScode\project\python\神经网络实验\MNIST\t10k-images.idx3-ubyte"
test_label = r"D:\VScode\project\python\神经网络实验\MNIST\t10k-labels.idx1-ubyte"

#input是数据，output是标签
input_train = (decode_idx3_ubyte(train_image)/256-0.5)*2
output_train = decode_idx1_ubyte(train_label)
input_test = (decode_idx3_ubyte(test_image)/256-0.5)*2
output_test = decode_idx1_ubyte(test_label)
print("数据集读取完成")



#in_test为shape[0]*784*1，shape[0]为图像的个数；input_test为shape[0]*28*28
#for i in range(input_train.shape[0]):
#    input_train[i] = input_train[i].reshape(1,-1)
in_train = input_train.reshape(input_train.shape[0],-1,1)

out_train, output_train = label_reshape(output_train)
#output_train是原始的标签，内容为5，8等；out_train是更改为10个神经元的标签，内容为0,0,1,0,0,0,0,0,0,0
'''
#output_train是原始的标签，内容为5，8等；out_train是更改为10个神经元的标签，内容为0,0,1,0,0,0,0,0,0,0
output_train = output_train.reshape(-1,1)
out_train = np.zeros((output_train.shape[0],10))
for i in range(output_train.shape[0]):
    if output_train[i] == 0:
        out_train[i][0] = 1
    elif output_train[i] == 1:
        out_train[i][1] = 1
    elif output_train[i] == 2:
        out_train[i][2] = 1
    elif output_train[i] == 3:
        out_train[i][3] = 1
    elif output_train[i] == 4:
        out_train[i][4] = 1
    elif output_train[i] == 5:
        out_train[i][5] = 1
    elif output_train[i] == 6:
        out_train[i][6] = 1
    elif output_train[i] == 7:
        out_train[i][7] = 1
    elif output_train[i] == 8:
        out_train[i][8] = 1
    elif output_train[i] == 9:
        out_train[i][9] = 1
'''


#in_test为shape[0]*784*1，shape[0]为图像的个数；input_test是原始的数据，为shape[0]*28*28
#for i in range(input_test.shape[0]):
#    input_test[i] = input_test[i].reshape(1,-1)
in_test = input_test.reshape(input_test.shape[0],-1,1)

out_test, output_test = label_reshape(output_test)
#output_test是原始的标签，内容为5，8等；out_test是更改为10个神经元的标签，内容为0,0,1,0,0,0,0,0,0,0
'''
#output_test是原始的标签，内容为5，8等；out_test是更改为10个神经元的标签，内容为0,0,1,0,0,0,0,0,0,0
output_test = output_test.reshape(-1,1)
out_test = np.zeros((output_test.shape[0],10))
for i in range(output_test.shape[0]):
    if output_test[i] == 0:
        out_test[i][0] = 1
    elif output_test[i] == 1:
        out_test[i][1] = 1
    elif output_test[i] == 2:
        out_test[i][2] = 1
    elif output_test[i] == 3:
        out_test[i][3] = 1
    elif output_test[i] == 4:
        out_test[i][4] = 1
    elif output_test[i] == 5:
        out_test[i][5] = 1
    elif output_test[i] == 6:
        out_test[i][6] = 1
    elif output_test[i] == 7:
        out_test[i][7] = 1
    elif output_test[i] == 8:
        out_test[i][8] = 1
    elif output_test[i] == 9:
        out_test[i][9] = 1
'''


start = time.perf_counter()                  #开始时计时

innum, midnum, outnum = (784,88,10)
w1 = np.array(np.random.rand(midnum, innum)-0.5)*2*0.01
w2 = np.array(np.random.rand(outnum, midnum)-0.5)*2*0.01
#w*0.01或*0.1是为了使初始权重尽可能波动小，经过实验，可以提高模型准确度（n=10，92%-->96%；n=1，82%-->90%）
b1 = np.array(np.random.rand(1, midnum)-0.5)*2
b2 = np.array(np.random.rand(1, outnum)-0.5)*2

'''----------------------------------------------模型训练----------------------------------------------'''
#np.seterr(divide='ignore',invalid='ignore')
print("正在训练...")

loss_figure = []                     #loss函数给出的原数据
loss_draw = []

#学习率
xite = 0.01
#训练循环数
n = 10
for j in range(n):
    for i in range(in_train.shape[0]):#输入训练集的数据量
        x = in_train[i].reshape(1, in_train.shape[1])
        z2, out, a2, a3 = forward(x, w1, w2, b1, b2)
        #w1, w2, b1, b2 = grad(out, output_train[1, :], a2, a3, w1, w2, b1, b2, z2, x, xite)
        w1, w2, b1, b2 = grad(out_train[i], out, a2, a3, w1, w2, b1, b2, z2, x, xite)
        #if i%10000==0:    print(z2)
        loss_figure.append(loss(out_train[i], out))
        if i==20000:    print(loss_figure[i+j*60000])#查看同一个样本在多次训练后的loss
    print((j+1)*100/n,'%')

end = time.perf_counter()              #结束时计时

dimension = 0                          #dimension表示维度，这里有10个维度，可以从0-9
loss_figure = np.array(loss_figure)
for i in range(loss_figure.shape[0]):
    loss_draw.append(loss_figure[i][0])
plt.scatter(range(len(loss_draw)), loss_draw, marker='o', s=2)    
#marker = 'o'表示用原点表示，还可以三角（5）等；s=2为点的大小

#原来result出来的train_num, 0, dimension，现在在loss中进行求和，因此不需要了
#train_num表示训练的次数，即loss图中的横坐标             
#dimension表示维度，这里有10个维度，可以从0-9
plt.show()
'''-------------------------------------------模型训练结束-----------------------------------------------'''


weight1 = w1
weight2 = w2
bias1 = b1
bias2 = b2


#softmax每次输出相同的结果时，1.数据集和w都要划到（-1，1），而不能（0，1）；2.out和out_hat的相减顺序要对
'''----------------------------------------------模型测试----------------------------------------------'''
right = 0
wrong = 0
hat = 0
for i in range(in_test.shape[0]):                   #测试
    x = in_test[i].reshape(1, in_test.shape[1])
    z2, out, a2, a3 = forward(x, w1, w2, b1, b2)
    for k in range(10):
        if out[0][k]==np.max(out):
            hat = k
            #print(hat)
            break
    if hat==output_test[i][0]:
        right = right + 1
    else:
        wrong = wrong + 1
    #if i==2000:                  #测试计算准确率模块是否正确
        #print('k=',k)
        #print('output_test[i][0]=',output_test[i][0])
        #print(k==output_test[i][0])
    
    '''
    if i==5000:                   #测试数组对应关系是否正确
        print(output_test[i][0])
        np.set_printoptions(suppress=True)        #设置print选项的参数
        print('out: ',out)
        plt.imshow(input_test[i], cmap='Greys_r')
        plt.show()
    '''
'''--------------------------------------------模型测试结束---------------------------------------------------'''


print(right,wrong)
print('正确率为:',right/(right+wrong)*100,'%')
print('training time =',int((end-start)/60),'min',int((end-start)%60),'s','=',end-start,'s')

